SELF-SIMILAR VOIDING SOLUTIONS OF A SINGLE LAYERED 
MODEL OF FOLDING ROCKS 
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Abstract. In this paper we derive an obstacle problem with a free boundary to describe the 
formation of voids at areas of intense geological folding. An elastic layer is forced by overburden 
pressure against a V-shaped rigid obstacle. Energy minimization leads to representation as a non- 
linear fourth-order ordinary differential equation, for which we prove their exists a unique solution. 
Drawing parallels with the Kuhn- Tucker theory, virtual work, and ideas of duality, we highlight the 
physical significance of this differential equation. Finally we show this equation scales to a single 
parametric group, revealing a scaling law connecting the size of the void with the pressure/stiffness 
ratio. This paper is seen as the first step towards a full multilayered model with the possibility of 
voiding. 

Key words. Geological folding, voiding, nonlinear bending, obstacle problem, free boundary, 
Kuhn- Tucker theorem 

AMS subject classifications. 34B15, 34B37, 37J55, 58K35, 70C20, 70H30, 74B20, 86A60 

1. Introduction. The bending and buckling of layers of rock under tectonic 
plate movement has played a significant part in the Earth's history, and remains of 
major interest to mineral exploration in the field. The resulting folds are strongly 
influenced by a subtle mix of geometrical restrictions, imposed by the need for layers 
to fit together, and mechanical constraints of bending stiffness, inter-layer friction and 
worked done against overburden pressure in voiding. An example of such a fold is seen 
in Figure [lTlj here the voiding is visible through the intrusion of softer material (dark 
in this figure) between the harder layers (shinier in the figure) which have separated 
while undergoing intense folding. 




Fig. 1.1. A photograph of a geological formation from Millock Haven, Cornwall, UK, demon- 
strating the formation of voids, visible by the intrusion of softer material, while harder layers undergo 
intense geological folding. Scale is approximately 5m across. 

Consider a system of rock layers, of constant thickness, initially lying parallel 
to each other that are then buckled by an external horizontal force, while being held 
together by an overburden pressure. If rock layers do not separate during the buckling 
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process it is then inevitable that sharp corners will develop. To see this, consider a 
single layer buckled into the shape of a parabola with further layers, of constant 
thickness, lying on top of this. Moving from the bottom layer upwards, geometrical 
constraints mean that the curvature of the individuals layer tightens until it becomes 
infinite, marking the presence of a swallowtail singularity [JJ. Beyond this singularity 
the layers interpenetrate in a non-physical manner. This process is illustrated in 
Figure |1.2[ showing how the layers would continue through the singularity if they 
were free to interpenetrate. 



Fig. 1.2. A close-up view of the propagation of a sine wave, demonstrating the physically- 
unrealisable swallowtail catastrophe. 

Models have dealt with these singularities by for instance limiting the number 
of layers [H [7] , using the concept of viscosity solutions [JJ , or postulating a simpli- 
fied geometry of straight limbs punctuated by sharp corners, as is observed in kink 
banding [TOUTS]. These approaches, however, disregard the resistance of the layers to 
bending, which is expected to be especially relevant close to the singularity. Here we 
therefore introduce the property of elastic stiffness into the modelling, and combine 
it with a condition of non-interpenetration. As a result, the layers will not fit together 
completely, but do work against overburdern pressure and create voids.The folding of 
rocks is a complex process with many interacting factors. In a multilayered model it 
is clear that work needs to be done to slide the layers over each other in the presence 
of friction, to bend the individual layers and finally to separate the layers (voiding). 
In order to understand the interaction between the process of bending and voiding we 
will not consider the effects of friction in this paper but will leave this to the subject 
of later work. 

The process of voiding is illustrated in Figure |1.3| which shows a laboratory 
recreation of folding rocks obtained by compressing laterally confined layers of paper. 
As we move through the sample, the curvature of the layers increase until a point 
is reached where the work against pressure in voiding balances the work in bending 
and the layers separate. A number of features of the voiding process can be seen in 
this figure. It is clear that the voids have a regular and repeatable form and that a 
typical void occurs when a smooth layer of paper separates from one which has a near 
corner-shape. 

In this paper we present a simplified energy-based model of voiding inspired by 



the processes observed in Fig. 1.3 The model consists of a single elastic layer with a 



vertical displacement w(x) forced downwards, and bent, into a corner-shaped obstacle 



of shape f(x) < w{x) by a uniform overburden pressure q (see Fig. 1.4 1. The corner 
is defined to have infinite curvature at the point, x = 0. For \x\ sufficiently large, the 
layer and obstacle are in contact so that w = f . However, close to x = the layer 
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and obstacle separate, leading to a single void for those values of x for which w > f. 
We study both the resulting shape of this elastic layer and the size of the voiding 
region. This investigation is the first part of a more general study of the periodic 



multi-layered voiding pattern seen in Fig. 1.3 



To study this situation we construct a potential energy functional V{w) for the 
system, derived in Section[2j which is given in terms of the vertical displacement w(x) 
and combines the energy Ub required to bend the elastic layer and the energy Uy 
required to separate adjacent layers and form voids. The potential energy function is 
then given by 

V = U B + Uy= - X t, 5/2 dx + q (w-f) dx, where w > f. (1.1) 

The resulting profile is then obtained by finding the minimiser of V over all suitably 
regular functions w > f. This constrained minimization problem is closely related to 
many other obstacle problems, as can be found in the study of fluids in porous media, 
optimal control problems, and the study of elasto-plasticity [5]. 

While obstacle problems are often cast as variational inequalities [12] , here we use 
the Kuhn- Tucker theorem for its suitability when interpreting the results physically. 
In Section [2] we prove various qualitative properties of constrained minimizers, and 
use the Euler-Lagrange equation to derive a fourth-order free-boundary problem that 
they satisfy. 

In addition, we show that stationarity implies that a certain quantity (the 'Hamil- 
tonian') is constant in any region of non-contact (Section [3]). This property extends 
the well-known property of constant Hamiltonian in spatially invariant variational 
problems, going back to Noether's theorem. However, we also give a specific inter- 
pretation of both the fourth-order differential equation and the Hamiltonian in terms 
of horizontal and vertical variations, with clear analogues with the concept of virtual 
work. Here horizontal and vertical variations define virtual displacements on the sys- 
tem, and the resulting ODEs describe the required load balance at a given point of 
a stationary solution. In Section |3.3| we show how integration of the Euler-Lagrange 
equation and the Hamiltonian gives vertical and horizontal force balances for the 
system, where individual terms can be identified with their physical counterpart. 




Fig. 1.3. A laboratory experiment of layers of paper constrained and loaded. In this figure 
the black lines are for illustrative purposes, and are produced by inserting a single black layer of 
paper between 25 layers of white. The resulting deformation shows the formation of voids when the 
imposed curvature becomes too high. Note the regular and repeatable nature of the voids. 
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Fig. 1.4. This figure shows the setup of the model discussed in this paper. An overburden 
pressure q forces an elastic layer w into another layer f with a corner singularity. 1+ and l_ define 
the first points of contact either side of the centre line. In this paper the layer is described both by 
Cartesian coordinates (x, w) measured from the centre of the singularity, and intrinsic coordinates 
characterised by arc length s and angle t/>. 



Section [4] gives a shooting argument that shows there exists a unique solution to 
this obstacle problem. These can be rescaled to form a one-parameter group, which 
gives the main result of Section [5} 

Theorem 1.1. Given k > so that f = k\x\, there exists a constant (3 — B(k) > 
such that for all q > and B > 0, the horizontal size of the void £ and the vertical 
shear force at the point of contact Bw xxx (£—) scales so that 

/on-1/3 8 /a\ 2 / 3 

£ = B(l) Bw xxx (£—) = —B- 



BJ xx " v ' (l + k 2 f/ 2 \B 

In Section[6]we show that these analytical results agree with the numerics, as well 
as with physical intuition. As the ratio of overburden pressure to bending stiffness 
becomes large, the size of the void tends to zero, giving a deformation with straighter 
limbs and sharper corners. By allowing the layers to form a void, the model is capable 
of producing both gently curving and sharp-cornered folds, without violating the 
elastic assumptions. Understanding this local behavior at areas of intensive folding 
may be seen as a first step to a multilayered model with the possibility of voiding. 

2. A voiding model close to a geometric singularity. 

2.1. The modelling. We consider an infinitely long thin elastic layer, of stiff- 
ness B, whose deformation is characterized by its vertical position w(x) as a function 
of the horizontal independent variable i£R. Overburden pressure, from the weight 
of overlying layers, acts perpendicularly to the layer with constant magnitude q per 
unit length. The layer is constrained to lie above the a V-shaped obstacle, defined by 
the function f(x) = k\x\, i.e. w > f. Although we appear to solve the problem for 
an infinitely thin layer, the analysis is the same for any layer of uniform thickness up 
to changes of stiffnes B. In all cases w(x) defines the ceterline of the layer, and f(x) 
defines the shape the layer would take in the absense of voids. This is only possible 
in this special case since / has straight limbs, and can therefore be propagated for- 
wards and backwards without change. The setup and parameters of the model are 



summarised in Fig. 1.4 
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The contact set of a function w is the set T(w) — {x € K : w(x) — /(x)}, 
the non-contact set T c (w) is its complement, and we define the two contact limits 
£ + = hxi{x > : u(x) = f(x)} and I- = sup{x < : u(x) = f(x)}. 

We now derive a total potential energy function for the system, described by the 
displacement w. 



2.1.1. Bending Energy. Classic bending theory (e.g. [T7J Ch. 1]) gives the 

B 
2' 



bending energy over a small segment of the beam ds as dUs = |k(s) 2 (1s, where n is 



curvature. Integrating over all s we find 

B r 2j B f°° w 2 xx ds B f°° w 2 xx , 
Ub = tt / n 2 ds=— j- ■ xx ^ —dx = — / - x ^rT7Kdx 



2 7-oo 2^(1 + ^)3^ 2^(1 + ^)5/2 

The quadratic dependence on w xx implies that a sharp corner has infinite bending 
energy. This is the basic reason why at any finite overburden pressure the elastic layer 
will show some degree of voiding. 

2.1.2. Work done against overburden pressure in voiding. The overbur- 
den pressure acting on the layer is q per unit length, therefore considering displace- 
ments w for which w > / the work done by overburden pressure in voiding is given 
by q(w — f) dx, and integrating over all x gives 

/oo 
(w - f)dx. 
-oo 

We see that if q is large, then Uy becomes a severe energy penalty. 

2.1.3. Total potential energy. The total potential energy function is the sum 
of bending energy and work done against overburden pressure, 

„ , X %^ dx + q / (w-f)dx (2.1) 



2 J_ x (1 + ^)5/2 



The solutions of the system are then minimizers of the energy functional (2.1 ) subject 
to the constraint w > f. 

A natural space on which to define V is the complicated-looking H 2 oc (R) fl (/ + 
L 1 (M)). Here H 2 oc (R) is the space of all functions with second derivatives in L 2 {K) 
for any compact set K C R. Finiteness of the first term in V requires (at least) 
w € H 2 oc (R), and well-definedness of the second term requires that w — f € L 1 (M). 
However, under the conditions w > f and V(w) < oo these conditions are automati- 
cally met, and therefore we will not insist on the space below. 

2.2. Constrained minimization of total potential energy. 

2.2.1. Properties and existence of minimizers. Before deriving necessary 



conditions on minimizers of (2.1) under the condition w > /, we first establish a few 
basic, but important, properties. These are that a constrained minimizer exists, is 
necessarily convex and symmetric, and has a single interval in which it is not in contact 
with the obstacle. We will prove uniqueness using different methods in Section |4j 

We write w& for the convex hull of w, i.e. the largest convex function v satisfying 
v < w. If w > f, then since / is convex, it follows that w# > f. 

Theorem 2.1. For any w, V(w^) < V(w), and any constrained minimizer w is 
convex. For all x G K, — k < w x (x) < k. 
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Proof. First we note that if w £ H 2 oc (R), also w* £ H 2 oc (R). Indeed, by consid- 
ering expressions of the form 



wf(x 2 )-w*(x 1 )= w* x (x)dx, 



it follows that the measure wf x is Lebesgue-absolutely continuous, and satisfies < 
w* x <\w xx \. Since w xx £ L 2 {R), it follows that w*. £ L 2 {R). Then w* £ H 2 (K) for 
all compact KcRby integration. 

Defining the set ft :— {x £ R : w#(x) — u(x)}, the function w# is twice differ- 
entiable almost everywhere on fi, with a second derivative w^ x equal to w xx almost 
everywhere on Q. On the c omp lement Q c , wf x = by [3 Theorem 2.1]. 



Substituting w# into (2.1 1 shows that V(w) > V(w^), with equality only if 
w& = w. Since w minimizes V, we have w = w&, and therefore w is convex. 

The restriction on the values of w x follows from the monotonicity of w x and the 
fact that w — / tends to zero at ±oo. □ 

As a direct consequence of Theorem |2.1| 

Theorem 2.2. The non-contact set T c (w) of a minimizer w is an interval con- 
taining x = 0, and for all x > £ + and x < £_ we have w(x) = f(x). 

Note that this statement still allows for the possibility that l± = ±oo. 

Proof. Suppose that x± 7 X2 > are such that w = f at x = x\ and at x = x-2.. By 
convexity of w we then have w — f on the interval [351,0:2]. If the contact set T(w) 
is bounded from above, then by the convexity of w, there exists e > and a > 
such that w(x) > a + (k + e)x for all x £ K, implying that Uyiw) = 00. Therefore 
T(w) D [0, 00) is an interval, and if it is non-empty, then it is necessarily extends to 
+00. Similarly, T(w) D (— 00, 0] is an interval, and if non-empty it extends to —00. 

Finally, note that x — can not be a contact point, since the condition w > f 
would imply that w £ ifj^R). Therefore the non-contact set T c (w) is an interval 
that includes x — 0. □ 

Theorem 2.3. Any minimizer w is symmetric, so that w(x) = w(—x). 

Proof. We proceed by using a cut-and-paste argument. If w is a minimizer, then 
it follows from Theorems |2.1| and |2.2| that w is convex and for all x > £+ and x < £_ , 
w(x) = f(x). Therefore w x (£±) = ±fc, and the intermediate value theorem states 
that there exists x £ (£-,£ + ) such w x (x) = 0. If V^^ ^iw) > V^ x ^(w), then we 
define the function 



w 



w(x + x) — k\x\ if x < 0; 
tu(— (x + x)) — k\x\ if x > 0. 



If Vt-oo^w) < Vji oo] (w), then we define w as 



w{— (x + x)) — k\x\ if x < 0; 
w(x + x) — k\x\ if x > 0. 



In either case w £ H 2 oc , w is symmetric, and V(w) < V(w). 

Since w is a minimizer, w solves a fourth-order differential equation in its non- 
contact set T(w) c (which includes x = 0; see (2.6) and Remark 2.9 1. By standard 
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uniqueness properties of ordinary differential equations (e.g. [5j), w and w are identical 
on both sides of x = 0, and remain such until they reach the constraint /. Therefore 
w = w and is therefore symmetric. □ 

Corollary 2.4. Since w is symmetric, £+ = —£- = £. 

Finally, these assembled properties allow us to prove the existence of minimizcrs: 

Theorem 2.5. There exists a minimizer ofV subject to the constraint w > /. 

Proof. Let w n be a minimizing sequence. By Theorems |2.1| and |2.3| we can 
assume that w n is convex and symmetric, and we therefore consider it defined on 
M + . By the convexity, since w n {x) — f{x) — > as x — > oo, the derivative w n , x 
converges to /'(oo) = k as x — > oo; therfore the range of w n , x is [0, k]. Since by 
convexity J °°(w n — /) dx > w n (0) 2 /2k, the boundedness of V{w n ) implies that w n {0) 
is bounded. 

From the upper bounds on w n x it follows that w n-xx is bounded in L (M. + ); 
combined with the bounds on iw n (0) and w n ^ x (0) = 0, this implies that a subsequence 
converges weakly in H 2 (K) to some w for all bounded sets K C [0, oo). Since therefore 
w n<x converges uniformly on bounded sets, it follows that w x (0) = and that 

w 2 r°° w 2 

lim inf / , n f x ^, dx> / ''. . .dr. 



— io (1 + <J 5/2 'Jo (l + ^) 5/2 

Similarly, uniform convergence on bounded sets of w n , together with positivity of 
w n — /, gives by Fatou's Lemma 



lim inf 

n— f oo 



/>oo POO 

' (w n - f)dx> (w- f) dx. 
o Jo 



Therefore V(w) < lim inf V(w n ), implying that w is a minimizer. □ 

2.2.2. The Euler- Lagrange equation. We now apply the Kuhn- Tucker theo- 
rem [HI pp. 249] to derive necessary conditions for minimizcrs of ( |2.1| su bject to the 
constraint w > /. Since any minimizer w is symmetric by Theorem |2.3[ we restrict 
ourselves to symmetric w, and therefore consider w defined on K + with the symmetry 
boundary condition w x (0) = 0. 

Theorem 2.6. Let q, B,k > 0. Define the set of admissible functions 

A={wef + H 2 (R+) n L 1 (R+) : w x (0) = 0} . (2.2) 



If w minimizes (2.1) in A subject to the constraint w > f, then it satisfies the sta- 
tionarity condition 

f 00 \ u w xx 5 w 2 xx w x 1 f°° 

for all (p £ H 2 (IR + ) (1 L 1 (M + ) satisfying tp x (Q) = 0, where n is a non-negative measure 
satisfying the complementarity condition J Q (w — f) dfi = 0. 

Proof. For the application of the Kuhn- Tucker theorem we briefly switch variables, 
and move to the linear space X := H 2 (R + ) n L 1 (M, + ), taking as norm the sum of the 
respective norms of H 2 and L . For any w € A, we define the void function v :— w—f, 
which is an element of X; the two constraints v x (0) := w x (0) — f x {0+) = —k and 
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v := w — f > are represented by the constraint Q(v) < 0, where Q : X — > Z 
Kxlx H 2 (R + ) is given by 

v x (0) + k 
G(v) := I -v x (0)-k 



We also define V(v) := V(v + /). 

If w satisfies the conditions of the Theorem, then the corresponding function 
v e X minimizes V subject to Q{v) < 0. The functionals V : X — »• M and 
are Gateaux diffcrentiable; since Q is affine, v is a regular point (see [T^l p. 248]) of 
the inequality Q{v) < 0. The Kuhn- Tucker theorem []3J p. 249] states that there 
exists a z* in the dual cone P* = {z* e Z* : (z* , z) > Vz e Z with z > 0} of the 
dual space Z* , such that the Lagrangian 

C(-):=V(-) + (G(-),z*) (2.4) 

is stationary at v; furthermore, (Q(v),z*) = 0. 

This stationarity property is equivalent to ( 2.3 ) . The derivative of V in a direction 



ip G X gives the left-hand side of (2.3); the right-hand side follows from the Riesz 
representation theorem |16l Th. 2.14]. This theorem gives two non- negative numbers 
Ai and A2 and a non-negative measure fx such ((0,6, it), z*) = \\a + A26 + J ud/j 
for all a, b £ R and m e I. Therefore (Q'(ip),z*) = —ftpd/jior any 93 g AT with 
^(0)=0. 

In addition, the complementarity condition (G(v),z*) = implies L vd/i = 
J™(w-f)d» = 0. □ 

This stationarity property allows us to prove the intuitive result that all minimiz- 
ers make contact with the support /: 

Corollary 2.7. Under the same conditions the non-contact set, T(w) c , is 
bounded, i.e. £ < 00. 

Proof. Assume that the contact set T(w) is empty, implying /i = 0. In ( |2.3[ ) 
take (p n ( x ) — n) for some <p € C^°{R) with J ip dx ^ 0. Since w — f € L 

and w I3 . € L 2 , we have 111^(0;) — >■ k as x — > 00; therefore, as n — > 00, the translated 
function 

n)xx 1 , \ 
(1 + w^.fl A 



converges weakly to zero in L , implying that the first term in (2.3), with ip = tp n , 



¥>n,xxdx= xx , {y + n)ip X x(y)dy 



(1 + ^)5/2-,- ^(1 + ^2)5/2 

vanishes in the limit n —¥ 00. The second term vanishes for a similar reason. In the 
limit n->cowe therefore find q J ip dx = 0, a contradiction. □ 

The boundedness of the non-contact set now allows us to apply a bootstrapping 
argument to improve the regularity of a minimizer w, and derive a corresponding 
free-boundary formulation: 
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Theorem 2.8. Under the same conditions as Theorem \2.b\ the function w has 
the regularity w G C°°(T(w) c ) (1 C 2 (R + ), w xxx is bounded, and w xxxx is a measure; 
the Lagrange multiplier [i is given by 



[i = q£5i + qH{ - - l)££, 



(2.5) 



where H is the Heaviside function, and _S? is one- dimensional Lebesgue measure. In 
addition, w and fi satisfy 



D 



10 



w 



35 w xx w x 



(l + u£)5/2 (l + W 2)7/2 2(1 + ^2)7/2 2(1 + ^2)9/2 



+ ? = M (2.6) 



Finally, w also satisfies the free-boundary problem consisting of equation (2.6 I on 
{Q,£) (with fi = 0), with fixed boundary conditions 



w x (0) = and w xxx (0) = 0, 
and a free-boundary condition at the free boundary x = t, 

w(£) — k£, w x {£) — k, and w xx (£) = 0. 



(2.7) 



(2.8) 



Before proving this theorem we remark that by integrating (2.6 1 we can obtain 
slightly simpler expressions. From integrating (2.6) directly, and applying (2.7), we 
find 



D 



(l+U)2)5/2 2 {l + W%yl 



' II 1 , + qx = qxH(x - i) for all x > 0. (2.9) 



By substituting the free boundary conditions at x = £ into (2.9) we also hnd that the 
limiting values of w xxx at x = £ are given by 



E ,(M = -(l + fc 2 ) 5/2 -k w xxx (£+)=0. 



B 



(2.10) 



In addition, by multiplying (2.9) by w xx and integrating we also obtain 

2~ (1 + w 2)5/2 + ^ XW * ~ w ) = ^WxxiO) ~qw(0). 



(2.11) 



Note that the right-hand side of (2.9) does not contribute to the the integral since 
Wxx = for all x > £. Substituting the boundary conditions (2.8), we derive the 
condition 



^Bw xx (0) 2 = qw(0), 



(2.12) 



so that the previous equation becomes 



B 



2 (1 + ^)5/2 



+ q(xw x — w) = 0. 



(2.13) 
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Proof of Theorem \2.8[ Once again we switch variables to the void function, 
v := w — f and define the functions 



9 = B 



(1 + (v x + fc) 2 ) 5 / 2 



and 



5„ vl x (v x + k) 
2 {1 + (v x + k) 2 )V 2 



by (|2.3|) we make the estimate 
gfxx = 



htp x + / Qj,-q)(p< ||/i||2||<Ps||2 + ||j"-g||Tv|MI 
<C(||vJ a || 2 + lbx||i), 
where the total variation norm ||^||tv is defined by 

Hrv := su P { / Cd" = C e CQR+), ||C|U < °°}- 



Setti ng = ^i, it follo ws th at g is weakly differentiable, and g x £ L 2 +L°°. From 

.9x|(£,oo) = and therefore £ I? . By 



Theorem 



2.2 



and Corollary 



2.7 



calculating g x explicitly, we may write 



(1 + (v x + fc) 2 )i g x + 5 



1 + (v x + fc) 2 



(2.14) 



GL 1 



Theorem 



2.1 



shows that (1 + (v x + fc) 2 ) 5 / 2 6 therefore e L 1 , so that 

Vxx € L°° , which in turn shows that v xxx £ L 2 by (2.14 1. We also see that since 



h x = - 



ZVxxVxxxjVx + fc) 7 v xx( v x + k) 2 

(1 + (v x + fc) 2 ) 7 / 2 (1 + + fc)2)7/2 + (1 + ( Vx + fe )2)9/2 > 



(2.15) 



ei 2 



we have ft.^ £ L 2 . We now look to similarly bound v xxxx . In the sense of distributions, 
we have 



g xx = -h x + fi — q, 



(2.16) 



and since h x has bounded support, this is an element of Ai, the set of measures with 
finite total variation. We can now write 



eh 1 



,3 



h-L( 4-1^5/2 I 5 3v xxx V xx(Vx + fc) + v xx 35 v xx (v x + fc) 2 

i 1 V j j ^ ^ 2 (1 + (v x + fc) 2 2 (1 + ( Wa + fc) 2 )3/2 



continuous and bounded 



Since u X xxa; has finite total variation, v xxx is bounded. Calculating (2.16) explicitly 
we find 



B 



10 



(v x + k)v xx v x 



(1 + (v x + fc) 2 ) 5 / 2 (1 + (tfc + fc) 2 ) 7 / 2 2 (1 + (v x + fc) 2 ) 7 /2 

35 v xx (v x + k) 2 



D 



2 (1 + (v x + fc) 2 ) 9 / 2 



+ q = I*- 
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Switching back to the orginal variable w = v + f gives (2.6| 



We now turn to (2.5). From the complementarity condition J(w — f)d/j, = 
we deduce that supp/i c T(w). Theorem 2.2 shows that w = / on [£,oo), and 
substituting this directly into (2.6) shows that /x|u,oo) — Q-^\ U,oo)- This proves that 
[l has the structure 

fi = aS e + qH( - - l)<£ 
some a > 0. To determine the value of a, take (p £ C°°([0, oo)) with bounded support, 



and such that ip = 1 in [0, £+ 1] . Then the weak Euler-Lagrange equation ( 2.3 ) reduces 
to a = q£. 

We now turn to the boundary conditions. The boundary condition w x (0) = 
is encoded in the function space, and the natural boundary condition w xxx (0) = 
follows by standard arguments. The conditions w{£) — k£, w x {£) — k, and w xx (£) = 
all follow from the contact at x = £. □ 

Remark 2.9. An identical argument gives that any minimizer on R, without assuming 



symmetry, satisfies the equation (2.6) on 



3. Hamiltonian, intrinsic representation, and physical interpretation. 

In this section we pull together an number of key results. First we calculate the 
Hamiltonian for the system and discuss its interpretation in a static setting. We then 
show that both the Hamiltonian and the Euler-Largrange equation for the system 
can be represented in a combination of cartesian and intrinsic coordinates, which 
allows us to intepret both objects physically. This physical interpretation shows the 
correspondence between the rigorous mathematical description of the problem, seen 
in Section [2] and a physical understanding of the system. 

3.1. The Hamiltonian. There is a long history of applying dynamical-systems 
theory to variational problems on an interval. Elliptic problems can thus be in- 
terpreted as Hamiltonian systems in the spatial variable x [15j . implying that the 
Hamiltonian is constant in space. 

We apply the same view here. The conserved quantity "H, which we again call 
the Hamiltonian, is obtained by considering stationary points of the Lagrangian L 



in (2.4) with respect to horizontal or 'inner' variations x 1— >• x + eip for some <p G H 



2 



This defines a perturbed function w E (x) :— w(x + e<p(x)), whose derivatives are 

w x (x) = w x (x + etp(x))(l + eip x {x)), 

w e xx (x) = w xx {x + eip(x))(l + e(p x (x)) 2 +w x (x + ecp(x))e(p xx (x). 

The requirement that the Lagrangian C is stationary with respect to such variations 
gives the condition 



B 



'JxxxxW x „W XXX W XX W X 5 W XX W X 35 w xx w x 



10- 



(1 + ^2)5/2 ( 1 + w 2)7/2 2(1 + ^2)7/2' 2(1 + ^)9/2 

+ (q - n)w x = 0. (3.1) 



Integrating this equation we find that the left-hand side of the expression 
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is constant in x, and the fact that it is zero follows from its value at x = and (2.12). 



By analogy to the remarks above we call the left-hand side above the Hamiltonian. 



Note that equation (3.1) is equal to (2.6) times w x . This is a well-known phe 



nomenon in Lagrangian and Hamiltonian systems, and can be understood by remark- 
ing that the perturbed function w e can be written to first order in e as w+eipw x ; there- 
fore this inner perturbation corresponds, to first order in e, to an additive ('outer') 
perturbation of ipw x . 



3.2. Intrinsic representation. Equations (2.9) and (3.2) can be written in 
terms of intrinsic coordinates, characterized by the arc length s, measured from the 
point of symmetry x = 0, and the tangent angle tp with the horizontal. First we recall 
the relevant relations between the two descriptions: 



4> s = K = W XX /(1 + W 2 X ) 3/2 , ds/dx = (1 + W 2 X ) 1/2 , 

cos if) — -j- — 1/(1 + wl) 1 ' 2 , smij: = ^ = w x /(l + w 2 x fl 2 . 



ds 



ds 



(3.3) 
(3.4) 



First we rewrite the integrated Euler-Lagrange equation, (|2.9|), as 
d 



B 



dx 





+ 







l B W xx 



2 + 



qx 



= qxH(x-£)(l 



,2\ 



and apply (|3.3|) and (|3.4|) to obtain 

"1 



b£- [A] 

dx 



Bif) s sin ip + qx — qxH(x — €) 



I ds n 
sec^y— = 0, 

dx 



which can also be written as 



BtJj ss cos ijj + 2^^s sin?/; + qx = qxH(x — £). 



Similarly, the integral (2.13) may be represented as 



1 



Bill 2 cos -0 + q(x tan tp — w) = 0. 



(3.5) 



(3.6) 



Following a similar process, the Hamiltonian (3.2) can be rearranged to 

,2 



(1 + w x ) dx 



l B W 



1 



2 (l + «£)3 (l + ufe) I 



.(1 + ^)1 

In intrinsic coordinates this gives the expression 

1 



- + qw = kqxH(x — €). 



Bi\) ss sin?/; Bip s cosip + qw = kqxH(x — €). 



(3.7) 



Note that equations (3.5) and (3.7) can be combined to give 



Bip ss + qx cos if) + qw sin if) = 2qxH(x — £) cosip. 



(3.8) 
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3.3. Physical interpretation in terms of force balance. The combination of 



3.1 



Cartesian and intrinsic coordinates that we have used allow us to identify terms of ( 3.5 1 
and (3.7 1 with their physical counterparts. Figure 

from s 



demonstrates the forces acting 
to s = s, together with a conveniently chosen 
Note that force balances are conveniently calculated for 



on a section of the beam 
area of pressurized matter 
the solid object consisting of the beam and the roughly triangular body of pressurized 
matter (indicated by the hatching) ; this setup allows us to identify the total horizontal 
and vertical pressure, exerted by q, as qx and q(w(x) — w(0)). 



q(w-w(0)) 



qw(0) 



BA 




Fig. 3.1. Left: forces on a section of the beam with pressurized matter. Right: small element. 



The small element of the beam shown in Fig. |3.1| indicates how the well-known 
relations from small-deflection beam theory between lateral load q, shear force F, and 
bending moment M, 



dF = qds and dM = Fds, 



extend into large deflections. We now use these expressions to identify the terms 



of (3.51 and (3.7) 



The equilibrium equation (2.6) was obtained by additively perturbing w, i.e. by 
replacing w by w + stp. This corresponds to a vertical displacement, which suggests 
that (2.6) can be interpreted as a balance of vertical load per unit of length. The 



integrated version (2.9) indeed describes a balance of the total vertical load on the 
rod from s = to s = s — i.e. the total vertical load on the setup in Fig. |3.1] — as we 
now show. 



We write equation (2.9) in the intrinsic- variable version (3.5) as 



shear force F 

(Btfj s ) s cos?/' - 
v v ' 

vertical component of F 




vertical component of P 



qx 

total vertical 
pressure 



qxH(x 



total vertical 
contact force 



Since by definition M = Btp s , the term Bip ss — (Bip s ) s is the normal shear force F, 
and the first term above is its vertical component. The term qx is the total vertical 
load exerted by the pressure q between s — and s — s (see Fig. 3.1 ), and qxH(x — t) 



is the vertical component of the contact force. Finally, the only remaining force with 
a non-zero vertical component is the axial force P at x, which can be interpreted as a 
Lagrange multiplier associated with the inextensibility of the beam. This suggests the 
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interpretation of the only remaining term in the equation as the vertical component 
of P, implying that we can identify P as 



(3.9) 



We can do a similar analysis of the Hamiltonian equation (3.1 ). Since this equa- 



tion has been obtained by perturbation in the horizontal direction, we expect that 
integration in space gives an equation of balance of horizontal load. In the same way 



we write the integrated equation (3.2) in intrinsic coordinates (see (3.7)) as 



(Bijj s ) s sinip 

" v ' 

horizontal 
component of F 



1 



Btp^cosip + q(w — w(0)) + qw(0) = kqxH{x — £) 



horizontal component of P 



total horizontal 
pressure 



horizontal load 
at s — 



total horizontal 
contact force 



Then we similarly can identify the first two terms as the horizontal components of 
the shear force and the axial load, while the last term is the horizontal component 
of the contact force. The remaining two terms are the horizontal component of the 
pressure q and the axial force at s = 0; the fact that this latter equals qw(0) is 
consistent with (3.9) when one ta kes (|2.12 1 into account. 

Note that the axial load P of (3.9), falling from a maximum compressive value at 



the centre of the beam to zero at the point of support, appears as a nonlinear term 
dependent on the bending stiffness B. Such terms are not normally expected in beam 
theory where, unlike for two-dimensional plates and shells, bending and axial energy 
terms are usually fully uncoupled. It comes about because of the re-orientation of the 
axial direction over large deflections. 

4. Existence, uniqueness, and stability of solutions of the Euler-Lagrange| 
equation. The Kuhn- Tucker theorem only provides a necessary condition for a min- 
imizer; it provides no information about existence of one or many solutions, or about 
the stability of a solution. We now develop a shooting argument that proves both exis- 
tence and uniqueness for the free-boundary problem (2.6 2.8). This shooting method 
also motivates a numerical algorithm in Section [6j 

Theorem 4.1. Given q > ; B > 0, and k > 0, there exists a function w and a 
scalar £ > that solve the free-boundary problem of Theorem^. 



Proof. We consider the ODE (3.8) as an initial value problem with ip(0) — and 
^s(O) = A, where w is coupled to %j) by (3.4) and w(0) = 0. Since minimizers w are 
convex (Theorem 2.1), we take A > 0. For small s > we have tp € (0, 7r/2), and 
therefore 



iJJs = A ■ 



< A 



^ S ec(i,( S '))x(s')+ 1 -(^')) 2 t^Hs') 



ds' 



n B 



I, sec 2 (ip(s'))x'dx' < A - \x 2 . 



2B 



Hence, for all A > there is a point at x 
therefore w xx (£) 



£(A) < y/2BA/q at which ip s = and 



0. From (3.6) we deduce that 



-B^s cos-0 + q(xtanip 



w) = -BA 2 - qw(0). 
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Therefore at x = I we have 



1 



BA 2 + q(w — w(0)) = qxw x , 



and since q(w — w(0)) > at x — I, and x — £ < yj2BA/q, it follows that 



q \2 



A 



3/2 



< W x (t). 



Now, consider the case of small A, so that w is also small. To leading order we then 
have 



so that 



w xxx + — x = 0, w(0) = w x (0) = 0, ^^(0) = A, 



= A —x 2 , and w x = Ax %:X 3 . 

2B ' 6B 



This implies that if A is sufficiently small, then 



w x = -Ay 



I2B 



-A < k, 



and conversely if A is sufficiently large, then 



w x {l) > 



T ( i 



B V 2 



A 



3/2 



> fc. 



Hence, by continuous dependence of the solution w on A, there is a value of A and a 
value of I for which 

UJ X (£) = k and w xx (£) = 0. 

If we now translate the function w by adding the constant k£ — w(£), then the resulting 
function w fulfills the assertion of the theorem. □ 

We now show that this solution is in fact unique. 

Theorem 4.2. The solution of the free-boundary problem of Theorem \2.8\ is 
unique. 

Proof. The proof uses a monotonicity argument. Let ip(x,A) be a solution 
of (3.8 ((written as a function of x) with ip s (0) — ip x (0) = A > 0. Let Ai < A 2 ; 
for small x, ip(x,Ai) < ip(x,A 2 ). Let 

x := sup{x > : tp(x, Ai) < tp(x, A 2 )} > 0. 

Since w(x) — w(0) = tan?/; it follows that 

w(x,A 1 )-w{0,A 1 )<w(x,A 2 )~w(0,A 2 ), for all < x < x. (4.1) 



Rewriting (2.11) in terms of tp x gives 



ds 



dx (1+Wa;) 5 / 2 y B cos 3 ip 



-BA 2 + q(w — w(0)) — qxt&nip 
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Using (4.1| we deduce that for all < x < x, i/j x (x,Ai) < ifj x (x, A 2 ), which implies 
that x — oo. 

Now assume that there exist two different solutions ip(x, Ai) and ^(a;, A 2 ), with 
associated points of contact £i and £2 such that i/K^ii^i) = ^(£2, A 2 ) = arctan/c. 
Since we have shown that A2) > ^(o;, Ai), it follows that £2 < £\ (see Fig. 4.1). 
Since < w x (i, Ai) < k for all < x < l\, we have 



w{£ 1 ,A l )-w{i 2 ,A l )<k{ti-i2). 



(4.2) 



Evaluating ( 3.6 1 at the free boundary for the solutions ^( • , A,) and the corresponding 




Fig. 4.1. The diagram shows the monotonicity argument used to prove uniqueness. If 
t\ , Ai) = V(^2i A2) = arctan k and ip(x, A2) > ip{x, Ai) for all x > 0, then I2 Kti- 



functions w; — w(-,Ai) gives 

q(wi{£i) - Wi (0) + ^BA^ = qk£ l: i = 1,2. 
Writing the difference as 

q[(w 2 (t 2 ) - w 2 (0)) - ( Wl (£ 2 ) - ^i(O))] + |(A| - A?) 



+ g [fc(ii - 4) - K(4) - wi(4))] - 0, 



we observe from (4.1 1 and (4.2) that the left-hand side is strictly positive. This 
contradicts the assumption of multiple solutions. □ 

5. Scaling Laws. We now see how the solutions of Section [4] can be written 
as a one-parameter group parameterized by q/B. Let £(q, B, k) be the length of the 
non-contact set T{w) of the solution w for that q, J5, and k, as defined in Section [4] 

Theorem 5.1. Given k > 0, there exists a constant j3 = fi{k) > such that for 
all q> and B > 0, 

-1/3 



l(q,B,k)=/3(±) . (5.1) 
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Proof. If we let x —: Xy, w —: Xu, and £ =: A/3, then the system (2.9) on (0, i) 
rescales to 



'mm 



(1+^)5/2 2(1+^)7/2 



" 7 "'"" +A 3 |y = on(0,/3). 



B- 



(5.2) 



By choosing A such that X 3 q/B = 1, the problem reduces to that of finding a w and j3 
such that 



M yyy 



2 

u yy u V 



{l + u 2jo/2-\ B {l + V ^y,2 +V = {) > u y (0)=0,u y ((3) = k, and %3/ (/3) = 0. 



Theorems |4 . 1 1 and 4.2 prove that for each k > there exists a unique pair (/?, u) that 
solve (5.2). Transforming back, (5.1) follows. □ 

Since «w(£) = -|(1 + fc 2 ) 5 / 2 £ (see ( |2.10[ >), it follows that 
Corollary 5.2. The shear force w xxx (£—) satisfies 

/3 /o\2/3 

w xxx [q,B,k](t-) = - 



(l + fc 2 )5/2 \B 



6. Numerical results. Here we provide some numerical results to support the 
analytical results seen in the previous section. The shooting method of the previous 
section suggests a numerical algorithm, by reducing the boundary value problem to 
an initial value problem, and shooting from the free boundary with the unknown 
parameter I. A one parameter search routine was made using MATLAB's built-in 
function fminsearch, which is an unconstrained nonlinear optimization package that 
relies on a modified version of the Nelder-Mead simplex method [T3] . 

Finding global minimizers in an unknown energy landscape can lead to an un- 



stable routine; however in this case the linearized version of (2.6) admits an analytic 



solution which provides a sufficiently accurate initial guess. Over the non-contact 



region equation (2.6) linearizes to 



w xxxx + — = 0. (6.1) 



Integrating (6.1) and applying the boundary conditions at the free boundary x = I 



gives the solution 



-2^+2^2 + ^(0), 



with the closed-form solution for 



I = 



3k B 



Figure |6.1| shows examples of solution profiles obtained in this manner. 

For fixed fc, the parameters £ and w xxx (£—) = —q£ can be calculated numerically 



for varying values of q/B, and the results are shown in Fig. 6.2 These numerical 



results agree with the behaviour expected. For fixed B, increasing q decreases the 
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-1 -0.5 0.5 1 



Fig. 6.1. Solution profiles for fixed q = 1, B = 1 and for increasing values of k 



a) b) 




Fig. 6.2. Numerical results supporting the scaling laws found for I and w xxx (£—) in Section^ 
results are shown for a fixed value of k = 0.75. (a) *'s show results found numerically for I against 
■g, the line plots f3 1//3 (b) *'s show results found numerically for w xxx (£—) against -J , the 

line plots - (1+fc /9 2)5/2 (f ) 2/3 . 
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size of the delamination, yet increases the vertical component of shear at delamination, 
Bw xxx (£—). Numerically fitting these curves to a function of the form /3(g) , we see 
that the results agree with the scaling laws found in the previous section, so that 

I = P (f )~ V3 , *>***{£-) = -(i + k 2 ) 5/2 P (|) 2/3 ■ 

Finally, Fig. |6.3| illustrates the dependence of j3 on k. 








Fig. 6.3. £ = fi{k) versus k. Here q = l and B = 1. 

7. Concluding remarks. The results of this paper show how elasticity, overbur- 
den pressure, and a V-shaped obstacle work together to produce one of the hallmarks 
of geological folds: straight limbs connected by smooth hinges. The model also gives 
insight into the relationship between material and loading parameters on one hand 
and the geometry and length scales of the resulting folds on the other. 

The model is of course highly simplified, and many modifications and general- 
izations can be envasiged. An important assumption is the pure elasticity of the 
material, and there are good reasons to consider other material properties of the 
layers. However, we believe the more interesting questions lie in other directions. 

One such question is role of friction between the layers, which was shown to be 
essential in other models of multilayer folding [TTJ [TOl 13 HH1 U\ ■ Since the normal stress 
between the layers changes over the course of an evolution, the introduction of friction 
will necessarily introduce history dependence, and the current energy-based approach 
will not apply. In this case other factors will also influence the behaviour, such as the 
length of the limbs, which determines the total force necessary for interlayer slip. 

An even more interesting direction consists in replacing the rigid, fixed, obstacle 
by a stack of other layers, i.e. by combining the multi-layer setup of pQ with the 
elasticity properties of this paper. A first experiment in that direction could be 
a stack of identically deformed elastic layers. An elementary geometric argument 
suggests that reduction of total void space could force such a stack in to a similar 



straight-limb, sharp-corner configuration, as illustrated in Fig. 7.1 This suggests that 
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Fig. 7.1. Sharp-angle, straight-limb folds give rise to fewer voids than rounded folds (after \1CX ). 

this phenomenon should also be visible in a stack of compressed layers, and we plan 
to consider this problem in future work E] . 
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